#联合尝试测试为20个点，少写无穷，成功
from gurobipy import *
#订单数据结构
# class Data:
#     timepiece = 10
#     multipleNum = 10
#     OrderNumber = 12 #假设订单0为depot，n+1为depot
#     orderNum = [0,1,2,3,4,5,6,7,8,9,10,11]  # 订单编号
#     customerNo = [0,1,2,3,4,5,6,7,8,9,10,0] #对应的顾客编号
#     #orderTime = []  # 订单到来时间
#     #prepareTime = []  # 需要的拣货时间
#     orderDemand = [0,10,10,20,20,20,20,10,10,30,10,0]  # 需求量
#     serviceTime = [0,10,10,10,10,10,10,10,10,10,10,0]  # 服务时间
#     dueTime = [1000,200,200,200,200,200,400,400,400,400,500,1000]    #最迟时间的软时间窗
#
# #客户数据结构
# class Customer:
#     allNumber = 11   #包含配送中心
#     # Cor_X_depot = 0
#     # Cor_Y_depot = 0
#     allCus = [0,1,2,3,4,5,6,7,8,9,10]
#     cor_X = [40,45,45,42,42,42,40,40,38,38,35]
#     cor_Y = [50,68,70,66,68,65,69,66,68,70,66]
#     distanceMatrix = [[]]
# class Vehicle:
#     vehicleNum = 3
#     vehicleCapacity = 60
class Data:
    timepiece = 10
    multipleNum = 15
    OrderNumber = 22 #假设订单0为depot，n+1为depot
    orderNum = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]  # 订单编号
    customerNo = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,0] #对应的顾客编号
    orderTime = [0,0,0,0,0,0,0,0,20,35,40,60,65,68,70,80,115,120,122,128,128,0]  # 订单到来时间
    prepareTime = [0,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,0]  # 需要的拣货时间
    orderDemand = [0,10,10,20,20,20,20,10,10,30,10,10,10,20,20,20,20,10,10,30,10,0]  # 需求量
    serviceTime = [0,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,0]  # 服务时间
    #dueTime = [1000,150,150,150,150,150,200,200,200,200,200,300,300,300,300,300,300,300,300,300,300,1000]    #最迟时间的软时间窗
    dueTime = [300, 150, 150, 150, 150, 150, 150, 150, 200, 200, 200, 200, 200, 200, 200, 300, 300, 300, 300, 300, 300,
               300]
    timeMatrix = [[]]  # 记录订单到达时间的区间 定义w(i_s)
    timeNum = []
#客户数据结构
class Customer:
    allNumber = 21   #包含配送中心
    # Cor_X_depot = 0
    # Cor_Y_depot = 0
    allCus = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
    cor_X = [40,45,45,42,42,42,40,40,38,38,35,35,25,22,22,20,20,18,15,15,30]
    cor_Y = [50,68,70,66,68,65,69,66,68,70,66,69,85,75,85,80,85,75,75,80,50]
    distanceMatrix = [[]]
#车辆数据结构
class Vehicle:
    vehicleNum = 6
    vehicleCapacity = 70
data = Data()
customer = Customer()
customer.distanceMatrix = [([0] * customer.allNumber) for p in range(customer.allNumber)]
for i in range(customer.allNumber):
    for j in range(customer.allNumber):
        temp = (customer.cor_X[i]-customer.cor_X[j]) ** 2 + (customer.cor_Y[i]-customer.cor_Y[j]) ** 2
        customer.distanceMatrix[i][j] = math.sqrt(temp)
        temp = 0
data.timeMatrix = [([0] * data.timepiece) for p in range(data.OrderNumber-2)]    #初始化，防止浅拷贝 timepiece列 OrderNumber行
for i in range(1,data.OrderNumber-1):
    for j in range(0, data.timepiece):
        if (data.orderTime[i] + data.prepareTime[i])/data.multipleNum < j :
            data.timeMatrix[i-1][j] = 1
data.timeNum = [0] * data.timepiece
for i in range(1, data.timepiece):
    for k in range(0,data.OrderNumber-2):
        data.timeNum[i] += data.timeMatrix[k][i]
vehicle = Vehicle()
BigM = 1e6
I = 5  #最少发车间隔的订单量
lam = 0.04  #发车次数下限比例
CPmixed = 100   #固定成本
ttload = 5   #装车时间
alpha = 1
beta = 1
gamma = 1
miu = 1
y = {}#发车时间
z = {}#发车时间订单对应
x = {}#路由
rsk = {}#车辆使用
tti = {}#服务客户时间
ttbsk = {}#返回时间
ttwi = {}#等待时间
# uusk = {}
model = Model('test')
#定义决策变量
#y_rs
for i in range(data.timepiece):
    for j in range(data.timepiece):
       if(i != j):
            name = 'y' + str(i) + '_' + str(j)
            y[i,j] = model.addVar(0
                                  ,1
                                  ,vtype = GRB.BINARY
                                  ,name = name)
#z_is
for i in range(data.OrderNumber):
    for j in range(data.timepiece):
        name = 'z' + str(i) + '_' + str(j)
        z[i,j] = model.addVar(0
                                  ,1
                                  ,vtype = GRB.BINARY
                                  ,name = name)
#x_sijk
for s in range(data.timepiece):
    for i in range(data.OrderNumber):
        for j in range(data.OrderNumber):
            for k in range(vehicle.vehicleNum):
                name = 'x' + str(s) + '_' + str(i) + '_' + str(j) + '_' + str(k)
                x[s,i,j,k] = model.addVar(0
                                  ,1
                                  ,vtype = GRB.BINARY
                                  ,name = name)
#r_sk
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        name = 'rsk' + str(s) + '_' + str(k)
        rsk[s,k] = model.addVar(0
                                  ,1
                                  ,vtype = GRB.BINARY
                                  ,name = name)
#tti
for i in range(data.OrderNumber):
    name = 'tti' + str(i)
    tti[i] = model.addVar(0
                          ,450
                          ,vtype=GRB.CONTINUOUS
                          ,name = name)
#ttbsk
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        name = 'ttbsk' + str(s) + '_' + str(k)
        ttbsk[s,k] = model.addVar(0
                          ,450
                          ,vtype=GRB.CONTINUOUS
                          ,name = name)
#ttwi
for i in range(data.OrderNumber):
    name = 'ttwi' + str(i)
    ttwi[i] = model.addVar(0
                          ,150
                          ,vtype=GRB.CONTINUOUS
                          ,name = name)
#uusk
# for s in range(data.timepiece):
#     for k in range(vehicle.vehicleNum):
#         name = 'uusk' + str(s) + '_' + str(k)
#         uusk[s,k]= model.addVar(0
#                        , 1
#                        , vtype=GRB.BINARY
#                        , name=name)
# for key in uusk.keys():
#     uusk[key].Start = 1
model.update()

#目标函数
obj = LinExpr(0)
for i in range(1,data.OrderNumber-1):
    for j in range(data.timepiece):
        obj.addTerms(alpha * j * data.multipleNum - alpha * data.orderTime[i] - alpha * data.prepareTime[i],z[i,j])


for s in range(data.timepiece):
    for i in range(data.OrderNumber):
        for j in range(data.OrderNumber):
            for k in range(vehicle.vehicleNum):
                if (i!=j):
                    obj.addTerms(beta * customer.distanceMatrix[data.customerNo[i]][data.customerNo[j]],x[s,i,j,k])

for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        obj.addTerms(gamma * CPmixed,rsk[s,k])

for i in range(1,data.OrderNumber-1):
    obj.addTerms(miu,ttwi[i])
model.setObjective(obj, GRB.MINIMIZE)


#相邻间隔最少订单数
for r in range(data.timepiece-1):
    for s in range(r+1,data.timepiece):
        model.addConstr(data.timeNum[s] - data.timeNum[r] >= I - BigM * (1 - y[r,s]),name='cons0')
#最少趟次约束
lhs = LinExpr(0)
for r in range(data.timepiece-1):
    for s in range(r+1,data.timepiece):
        lhs.addTerms(1 , y[r,s])
model.addConstr(lhs >= data.timeNum[data.timepiece-1] * lam,name = 'cons2')
#订单i在s时刻发出的必要条件是该时刻拣选完成
for i in range(1,data.OrderNumber-1):
    for j in range(data.timepiece):
        model.addConstr(z[i,j] <= data.timeMatrix[i-1][j],name = str(i)+'_'+str(j)+'biyao')
#订单i需要被配送
for i in range(1,data.OrderNumber-1):
    lhs = LinExpr(0)
    for j in range(data.timepiece):
        lhs.addTerms(1,z[i,j])
    model.addConstr(lhs == 1,name = 'onetime'+'_'+str(i))
#z和y之间的关系1
for r in range(data.timepiece-1):
    for s in range(r+1,data.timepiece):
        lhs = LinExpr(0)
        for i in range(1,data.OrderNumber-1):
            lhs.addTerms(1,z[i,s])
        model.addGenConstrIndicator(y[r,s] , 1 , lhs>=1 , name = 'relationship1'+str(r)+'_'+str(s))

#z和y之间的关系2
for i in range(1,data.OrderNumber-1):
    for s in range(1,data.timepiece):
        lhs = LinExpr(0)
        for r in range(0,s):
            lhs.addTerms(1,y[r,s])
        model.addGenConstrIndicator(z[i,s] , 1 , lhs==1,name = 'relationship2'+str(i)+'_'+str(s))
#起始区间流平衡
lhs = LinExpr(0)
for j in range(1,data.timepiece):
    lhs.addTerms(1,y[0,j])
model.addConstr(lhs == 1, name='flow_conservation_0')
#中间时间段的流平衡

for h in range(1,data.timepiece-1):
    expr1 = LinExpr(0)
    expr2 = LinExpr(0)
    for i in range(h-1):
        expr1.addTerms(1,y[i,h])
    for j in range(h+1,data.timepiece):
        expr2.addTerms(1,y[h,j])
    model.addConstr(expr1==expr2,name='flow_conservation_'+str(h))
    expr1.clear()
    expr2.clear()
#终止区间流平衡
lhs = LinExpr(0)
for j in range(0,data.timepiece-1):
    lhs.addTerms(1,y[j,data.timepiece-1])
model.addConstr(lhs == 1, name='flow_conservation_last')

#逻辑约束,需求被一辆车服务一次
for s in range(data.timepiece):
    for i in range(1,data.OrderNumber-1):
        lhs = LinExpr(0)
        for k in range(vehicle.vehicleNum):
            for j in range(1,data.OrderNumber):
                if (i != j) :
                    lhs.addTerms(1, x[s,i,j,k])
        model.addConstr(lhs == z[i,s],name = 'luoji1' + str(i) + '_' + str(s))
        lhs.clear()
# for s in range(data.timepiece):
#     for j in range(1,data.OrderNumber-1):
#         lhs = LinExpr(0)
#         for k in range(vehicle.vehicleNum):
#             for i in range(0, data.OrderNumber-1):
#                 if(i != j):
#                     lhs.addTerms(1,x[s,i,j,k])
#         model.addConstr(lhs == z[j,s],name = 'luoji2' + str(j) + '_' + str(s))
#         lhs.clear()
#只有当s时刻有订单发出时，s时刻才有车辆被使用
for s in range(data.timepiece):
    lhs = LinExpr(0)
    for i in range(1,data.OrderNumber-1):
        lhs.addTerms(1, z[i,s])
    for k in range(vehicle.vehicleNum):
        model.addConstr(lhs >= rsk[s,k],name = 'vehicle' + str(s) + '_' +str(k))
    lhs.clear()
#起点流平衡
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        lhs = LinExpr(0)
        for j in range(1,data.OrderNumber):
            if (j != 0):
                lhs.addTerms(1,x[s,0,j,k])
        model.addConstr(lhs == rsk[s,k],name = 'beginflow' + str(s) + '_' + str(k))
        lhs.clear()
#中间点流平衡
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        for h in range(1,data.OrderNumber - 1):
            expr1 = LinExpr(0)
            expr2 = LinExpr(0)
            for i in range(data.OrderNumber-1):
                if (h != i):
                    expr1.addTerms(1,x[s,i,h,k])
            for j in range(1,data.OrderNumber):
                if (h != j):
                    expr2.addTerms(1,x[s,h,j,k])
            model.addConstr(expr1 == expr2,name = 'middleflow' + str(s) + '_' + str(k) + '_' + str(h))
            expr1.clear()
            expr2.clear()
#终点流平衡
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        lhs = LinExpr(0)
        for i in range(data.OrderNumber - 1):
            if (i != 0):
                lhs.addTerms(1,x[s,i,data.OrderNumber - 1,k])
        model.addConstr(lhs == rsk[s, k], name='endflow' + str(s) + '_' + str(k))
        lhs.clear()
# #x_sijk,r_sk的关系
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        lhs = LinExpr(0)
        for i in range(data.OrderNumber):
            for j in range(data.OrderNumber):
                if(i != j):
                    lhs.addTerms(1,x[s,i,j,k])
        model.addGenConstrIndicator(rsk[s,k], 1, lhs >= 1, name='relationshipxr1' + str(s) + '_' + str(k))
        model.addGenConstrIndicator(rsk[s, k], 0, lhs == 0, name='relationshipxr2' + str(s) + '_' + str(k))
        lhs.clear()
#载重约束
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        lhs = LinExpr(0)
        for i in range(1,data.OrderNumber-1):
            for j in range(data.OrderNumber):
                if(i != j):
                    lhs.addTerms(data.orderDemand[i],x[s,i,j,k])
        model.addConstr(lhs <= vehicle.vehicleCapacity, name = 'capacity_vehicle' + str(s) + '_' + str(k))
        lhs.clear()
#每趟次第一个客户的到达时间约束
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        for j in range(1,data.OrderNumber):
            model.addConstr(customer.distanceMatrix[data.customerNo[0]][data.customerNo[j]]+ s * data.multipleNum + ttload - BigM * (
                    1 - x[s,0,j,k]) - tti[j] <= 0, name = 'time' + str(j))
#中间客户的到达时间约束
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        for i in range(data.OrderNumber-1):
            for j in range(1,data.OrderNumber):
                if (i != j):
                    model.addConstr(customer.distanceMatrix[data.customerNo[i]][
                                        data.customerNo[j]] + tti[i] + data.serviceTime[i] - BigM * (
                                                1 - x[s, i, j, k]) - tti[j] <= 0, name='time' + str(j))
#每趟次返回配送中心的时间约束
for s in range(data.timepiece):
    for k in range(vehicle.vehicleNum):
        for i in range(data.OrderNumber-1):
            model.addConstr(customer.distanceMatrix[data.customerNo[i]][
                                data.customerNo[data.OrderNumber-1]]+ tti[i] + data.serviceTime[i] - BigM * (
                                    1 - x[s, i, data.OrderNumber-1 , k]) - ttbsk[s,k] <= 0, name='time' + str(data.OrderNumber-1))
#时间窗惩罚
for i in range(1,data.OrderNumber-1):
    model.addConstr(ttwi[i] - tti[i] + data.dueTime[i] >= 0,name = 'timewindowviolate'+ '_' + str(i))
#返回时间超过准备出发时间时，车辆不可用
for k in range(vehicle.vehicleNum):
    for t in range(data.timepiece-1):
        for s in range(t+1,data.timepiece):
            #model.addConstr(ttbsk[s,k]-t * data.multipleNum - BigM * rsk[s,k] >= 0,name = 'useable1' + '_' + str(s) + '_' +str(k))
            model.addConstr(ttbsk[t,k] -s * data.multipleNum - BigM + BigM * rsk[s,k] <= 0, name = 'useable2' + '_' + str(s) + '_' +str(k))
# 导出模型
model.write('Test.lp')
# 求解
model.optimize()
# 打印结果
print("\n\n-----optimal value-----")
print(model.ObjVal)

for key in y.keys():
    if (y[key].x>0):
        print(y[key].VarName + ' = ', y[key].x)
for key in z.keys():
    if (z[key].x>0):
        print(z[key].VarName + ' = ', z[key].x)
for key in x.keys():
    if (x[key].x>0):
        print(x[key].VarName + '=', x[key].x)
for key in ttbsk.keys():
    if (ttbsk[key].x>0):
        print(ttbsk[key].VarName + '=', ttbsk[key].x)